# -*- coding: utf-8 -*-
"""
Created on Sat Nov 14 15:16:48 2020

@author: 310
"""
import numpy as np

def func(x):
    if x==0:
        return 1
    return np.sin(x)/x

def Romberg(a,b,e,func):
    h=b-a
    k=1
    T1=h/2*(func(a)+func(b))
    
    while True:
        S=0
        x=a+h/2
        while x<b:
            S=S+func(x)
            x=x+h
            
            
        T2=T1/2+h/2*S
        S2=T2+(T2-T1)/3
        
        if k==1:
            k=k+1
            h=h/2
            T1=T2
            S1=S2
            continue
        else:
            C2=S2+1/15*(S2-S1)
            if k==2:
                C1=C2
                k=k+1
                h=h/2
                T1=T2
                S1=S2
                continue
            else:
                R2=C2+1/63*(C2-C1)
                if k==3:
                    R1=R2
                    C1=C2
                    k=k+1
                    h=h/2
                    T1=T2
                    S1=S2
                    continue
                else:
                    if abs(R2-R1)<e:
                        return R2
                    else:
                        R1=R2
                        C1=C2
                        k=k+1
                        h=h/2
                        T1=T2
                        S1=S2
            
            
a=0
b=10
e=1e-10
result=Romberg(a,b,e,func)
from scipy import integrate
R=integrate.quad(func,a,b)[0]
print("err:%s"%(abs(R-result)))

                        